Tailoring THz radiation by controlling tunnel photoionization events in gases 
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Applications ranging from nonlinear terahertz spectroscopy to remote sensing require broadband 
and intense THz radiation which can be generated by focusing two-color laser pulses into a gas. 
In this setup, THz radiation originates from the buildup of the electron density in sharp steps of 
attosecond duration due to tunnel ionization, and subsequent acceleration of free electrons in the 
laser field. We show that the spectral shape of the THz pulses generated by this mechanism is 
determined by superposition of contributions from individual ionization events. This provides a 
straightforward analogy with linear diffraction theory, where the ionization events play the role of 
slits in a grating. This analogy offers simple explanations for recent experimental observations and 
opens new avenues for THz pulse shaping based on temporal control of the ionization events. We 
illustrate this novel technique by tailoring the spectral width and position of the resulting radiation 
using multi-color pump pulses. 



INTRODUCTION 



Recently, THz technology has attracted significant 
attention, because it can provide unique analytical 
and imaging tools in nonlinear and time-domain spec- 
troscopy, remote sensing, biology and medicine, security 
screening, information and communication technology 
(see e.g.[l|-lj|). Among the various THz sources, enrploy- 
ing two-color ionizing femtosecond pulses in gases |5l-[l4 
stands out through absent damage threshold, low phonon 
absorption and zero interface reflection. The two-color 
scheme provides striking performance in terms of band- 
width that can exceed 100 THz and thus even covers far- 
and mid-infrared, as well as high peak fields up to the 
MV/cm range. 

In the standard realization of this method a femtosec- 
ond laser pulse and its second harmonic are collinearly 
focused into a gas, such as ambient air, yielding power- 
ful THz radiation with well-defined electric field values. 
Initially, this process has been explained by four-wave 
mixing rectification via third-order nonlinearity [5|, Il2| . 
However, a threshold for THz gene ration [6|,ll0j as well as 
plasma current measurements J13l | indicated that plasma 
formation plays an important role in the generation pro- 
cess. Besides, the Kerr nonlinearity is far too small 
[71, |l5| to explain the observed THz field strength in 
such a setup. In order to solve this discrepancy, Kim et 
al. [ll|, UM have put forward the so-called photocurrent 
mechanism. In this mechanism tunneling ionization and 
subsequent electron motion produce a quasi-DC photoin- 
duced current, which in turn emits THz radiation. The 
quasi-DC plasma current can be efficiently produced in a 

H4l |. in multi- 
or few-cycle 



Results pubhshed so far [13, E, E H [SHl provide 
several important features of the THz spectrum. How- 
ever, a general framework for analysis, control and design 
of THz radiation is still missing. Similar to pulse shaping 
techniques in the optical range, the broadband control of 
THz waveforms and its spectral properties is highly desir- 
able for many applications, while the existing techniques 
are applicable only for a narrow bandwidth of few THz 
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gas irradiated by two-color laser pulses 
color pump schemes [l6| and in chirped 

pulses [il,[i3|. 



In this paper we advance the understanding of the fun- 
damental physics in THz generation in gases as well as 
examine its basic mechanism associated with plasma for- 
mation. We demonstrate that THz generation in gases is 
intrinsically connected to optically-induced stepwise in- 
crease of the plasma density due to tunneling ionization. 
New spectral components are emitted in a discrete set of 
attosecond-scale, ultra-broadband bursts associated with 
these ionization events; in the spectral representation 
THz emission results from an interference of these ra- 
diation bursts. 

Using appropriate field shapes of the pump pulses, e.g. 
tuning the frequencies and phases of two- or multi-color 
pulses, the temporal positions of the ionization events can 
be controlled. In turn, this allows to design the emitted 
radiation within a broad spectral range of more than 100 
THz. This new concept of prescribing the phase and am- 
plitude of each elementary contribution responsible for 
THz emission has a remarkable analogy in linear diffrac- 
tion theory: The produced spectrum can be understood 
like the far-field diffraction pattern of a grating, where 
the ionization events play the role of the slits. A compre- 
hensive (3-|-l)-dimensional non-envelope model of propa- 
gation is used to verify by means of numerical simulations 
that the simple picture drawn above captures the essen- 
tial physics of THz generation. 



STRUCTURE OF THE THZ SPECTRUM 

We start our analysis in the so-called local current 
(LC) limit, that is, we consider a small volume of gas 
irradiated by the strong pump field E{t) [see Fig. [TJa)]. 
Electrons produced due to ionization build up a current 
J{t) which, in turn, produces an electromagnetic wave 
E'-'{t) c>c dtJit) [27|. We focus here on the pump in- 
tensities where tunnel ionization is the dominant ioniza- 
tion process. In this regime, the free electron density 
p[t) increases stepwise in short attosecond-scale ioniza- 
tion events corresponding to maximum amplitudes of the 
pump field at times t„ [see Fig. [TJc)]. This stepwise in- 
crease of the free electron density p{t) was recently con- 
firmed experimentally [23, [29J . The ionization events are 
well separated from each other and thus p{t) can be writ- 
ten as a sum over contributions from all ionization events: 
p{t) = ^„SpnHn{t). Each ionization event has a well 
defined amplitude Spn and shape Hn{t). Since all ioniza- 
tion events have almost the same shape [see green filled 
curve in Fig. [UJc)], one can assume il„(i) ~ H{t — tn), 
where iJ is a "smoothed" step function. An approximate 
analytical formula for H{t) is given in the Methods. 

In Fig. [TJb-e) an exemplary two-color pump field 
E{t) — £ {t)[cos (ujQt) + r cos (2a;ot + 6*)] is considered, 
where ujq = 2tt x 375 ps^^ is the fundamental frequency 
(corresponding to the wavelength of 800 nni) , r — \/0.2 
and 6 — 7r/2 are relative phase and amplitude of the 
second harmonic (SH) component. For such a two-color 
pulse (and assuming r < 1) the positions of ionization 
events are given by wo^n ~ tto — 2(— l)"r sin0 [15|. f (t) 
is a Gaussian pulse envelope shape with full-width at 
half-maximum (FWHM) pulse duration t^ — 40 fs. 

With the assumption that electrons are born with zero 
initial velocity, the plasma current J(i) [blue curves in 
Fig. [Hb,c)] is governed by 



dtJ{t)^-lJ{t)^^E{t)p{t), 
m 



(1) 



where 7 is the electron-ion collision rate (7 = 5 ps~^ 
at atmospheric pressure) and m, q are electron mass and 
charge [ll|, ll5| . Then, in Fourier domain the generated 
electromagnetic wave in the LC limit E'^{uj) ex ujJiuj) 
reads (see Methods) 



E-'iu) 



E 



[A„(^)-C„B(c^)]e*' 



(2) 



Here, A„(cj) = qgujSpnT[H{t)vf{t + t„)], C„ = 
qSpnVf{tn), and B{uj) = gujF[H{t)e^'^^], with F denot- 
ing the Fourier transform, Vf{t) — — J_ E{T)e^^'^^^'dT 
is the free electron velocity, and 17 is a constant. 

Eq. ^ allows us to identify the impact of each ion- 
ization event. The interpretation of Eq. ^ is straight- 
forward when we consider the spectral dependence of the 



amplitudes An{uj) and -B(w), which can be calculated an- 
alytically in a reasonable approximation (see Methods): 



B{uj) ex 



w — 17 



■ exp 



(3) 



where cr = 10 fs^^ is the spectral width of a typical ion- 
ization peak [see thin black curves in Fig. [ijc)]. In con- 
trast to B{u]), the amplitudes A„(a;) depend nontrivially 
on the pump pulse. Examples of A„(a;) and B{ui) are 
shown in Fig. [ljd,e). i?(w) features two clearly sepa- 
rated frequency scales determined by 7 and a. On the 
short frequency scale w < 7 (up to 5 meV or ^ — 2 THz 
for atmospheric pressure), B{ui) ex uj, while for larger fre- 
quencies B{uj) depends only slowly on lu (on the scale of 
- 10 eV). 

If we consider the range from a few to approximately 
100 THz, Eq. ([2]) can be significantly simplified. In this 
range Biuj) can be considered as a constant, and addi- 
tionally the An{u}) are negligible as seen from Fig. 1(d). 



E-^ {<^) ^ B{uj)Y,Cne 



^iujtt 



OC ^ Cn^ 



^iujtji 



(4) 



In other words, the resulting THz spectrum is a simple 
linear superposition of the contributions _E;[ = C„e*"*" 
and its structure is mainly governed by interference due 
to the spectral phases e*"*". This interference is illus- 
trated in Fig. [TJe) [cf. also Fig. [2Kb)]. Note that in tem- 
poral domain Eq. Q corresponds to a set of sharp peaks 
localized near i„, each of them having the amplitude C„ 
and the shape dt [H{t - i„)e"''**"*"'] ■ Thus, Eq. Q de- 
scribes sharp attosecond-long bursts of radiation local- 
ized near the ionization events. 

Further on, the same simple superposition principle 
is also valid for higher frequencies in the more general 
Eq. ([2]). Indeed, Fig. [T^d) shows that constructive spec- 
tral interference in Eq. ([2|) appears also for frequencies 
UJ ~ /wo, where I is integer. This finding is consistent with 
recent experimental results [29| , where optical harmonics 
generated due to tunneling-ionization-induced modula- 
tions of the electron density were observed. Therefore, 
our superposition principle may also offer a complemen- 
tary understanding of Brunei radiation 29|, [30| . 

In fact, Eq. Q shows that THz generation can be un- 
derstood completely analogous to far-field interference 
patterns produced by a diffraction grating 31|, where 
the ionization events play the role of the slits. Each 
of such "slits" produces a broad "secondary wave" , and 
their interference results in narrow lines, according to the 
Huygens-Fresnel principle 3]|. THz radiation then cor- 
responds to the zeroth diffraction order. In contrast to 
diffractive gratings, interference here takes place in the 
frequency domain and not in position space. 

In a somewhat oversimplified picture, one can explain 
the spectrum shown in Fig.[IJd) analogous to interference 
patterns produced by diffraction gratings of equidistant 
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FIG. 1. Mechanism and spectral properties of terahertz generation by ionizing two-color pulses, (a) Schematic setup, (b) 
two-color pump field E{t) (red curve) and current J{t) (blue curve) induced by tunneling-ionization, (c) step-wise modulation 
of the free electron density (filled green curve) and ionization rate W{t) (thin dotted black curve), besides, few cycles of the 



electric field E(t) (red curve) and the current J(f) (blue curve) are shown, (d) Spectrum \E 



vs. frequency i^ = a;/27r on a 



large frequency scale (red curve). Each ionization event at time i„ induces a broadband radiation described by the form- factors 
B{v) (black curve) and An{i^) (blue curve). The sum of all partial contributions in Eq. ((Sjl yields the resulting spectrum 
E'' {v). (e) Partial contributions HeEn{v) from the nth ionization event for —16 < n < 16 (grey curves) and the total spectrum 
YleE^iv) (filled red curve) in the THz to mid-infrared spectral region {IraE'^ii') is negligible here). The black curve shows 
B{y). The decay of E'' (u) for u >20 THz is solely due to destructive interference of different E:^{u). In contrast, in the shaded 
region {u < 5 THz) E'^ (u) decays because of decrease of B{v). (f) Spectrum \E'\v)\'^ vs. input pulse duration rin. As the 
number of ionization events decreases with Tin, the spectral width increases, as also confirmed by 3D simulations for rin = 20, 
40, and 80 fs (black curves). 



slits, E'^{uj) ^ sm{Nuj5t/2)/ sm{u}St), where N is the 
number of ionization events and St ~ tt/wq is the time 
interval between two subsequent ionization events (see 
Auxiliary Material). The width of the spectral line is 
inversely proportional to N, that is, to the input pulse 
duration [see Fig. [Ijf)]. Such inverse dependence of the 
THz spectral width on the pump duration was already 
observed experimentally [7|. It is also compatible [black 
lines in Fig. [Uf)] with comprehensive 3D simulations, 
where propagation and transverse effects are accounted 
for (see Methods). 



CONTROL OF THE THZ SPECTRUM 

Understanding the spectrum of the radiation emitted 
by the plasma current J as an interference pattern pro- 
vides a deeper insight into the mechanism responsible for 
the THz generation. Beyond that, this understanding 
provides us with the possibility to shape the emitted ra- 
diation by suitable engineering of the pump field. Here, 
we will restrict ourself to the THz and far-infrared do- 



main described by Eq. ^. The latter indicates that, in 
order to control the THz spectral shape, we have to tar- 
get the times of the ionization events i„ and the values 

of C„ = qSpnVf{tn). 

Figure [5] shows several examples of THz generation 
from simple variations of the two-color scheme. Varia- 
tions of the relative phase between fundamental and 
second harmonic component of the pump pulse have a 
strong influence on the THz spectral shape [see cases (1) 
and (2)], in addition to the known dependency of the 
THz yield [l3J . This observation can be readily explained 
within our linear interference framework. For 9 = (case 
1) the coefficients C„ are not sign definite [red and blue 
circles in Fig. [2la)]. Thus, for w = different contri- 
butions in Eq. ^ cancel each other, that is, destructive 
interference dominates. Nevertheless, with increasing fre- 
quency id the destructive interference is (partially) com- 
pensated by the spectral phases e"^*". This behavior is 
shown in Fig.[2ljb), where contributions from the ioniza- 
tion events Ejl are visualized in the complex plane for 
different frequencies w. Consequently, the spectral max- 
imum is shifted from zero [see Fig. [5Jc)]. In contrast, for 
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FIG. 2. Interference of contributions from different ionization events (a,b), THz spectrum (c) and temporal shape (d) for 
different values of the phase 9 between the fundamental and second harmonics (cases 1, 2) as well as of the frequency shift 
Sv of the second harmonic (cases 3, 4; Si/ = =p20 THz, = 7r/4). (a) Times of ionization events t„ (center of circles) and 
amplitudes C„ [radii of circles, red (blue) color indicates positive (negative) sign]. The gray curves show the electric field 
amplitude E'^{t) of the pump pulse, (b) Values of the summands in Eq. Q E:^{u) in the complex plane for u = (red circles) 
and 1/ = 25 THz (green circles). For the cases (1) and (3) destructive interference occurs at ;/ = and is partially compensated 
at larger frequencies. (c),(d) Corresponding THz spectral shapes and temporal profiles obtained from Eq. ((4]) (X^n^"' filled 
green curve) and from 3D simulations at the linear focus z = 5 mm (dashed yellow curves). Note that the amplitudes for case 
(1) are about one order of magnitude smaller in comparison to the other cases. 



6 = tt/2 (case 2) all coefficients C„ are positive. This 
obviously ensures that at w = constructive interference 
dominates, and the phase factor e'"'" can lead only to 
destructive interference when w is increased. Therefore, 
the spectral maximum is located at zero frequency in 
this case. Moreover, the coefficients C„ are much larger 
in amplitude, which explains the increased THz yield for 
9 = tt/2. 

We verified this behavior by comprehensive (3+1)- 
dimensional (3D) numerical modeling (see Methods), 
which was recently successfully compared to experi- 
ments [l^. We consider a Gaussian input beam with 
width wq = 100 ^m and duration (FWHM) Tin = 40 fs 
for the 800 nm fundamental pump pulse (pulse energy 
~ 300 /LtJ). The energy of the second harmonic at 400 nm 
is chosen as 7.5 % of the fundamental, with duration and 
width smaller than the values for the fundamental by a 
factor ^/2 (similar to typical experimental conditions). 
These pulsed input beams are focused with / = 5 mm 
into argon at atmospheric pressure. In contrast to the 
simple LC model, the phase angle 9 between the two 
components shifts during propagation. For our particu- 
lar configuration, an initial value ^in = 7r/4 guarantees 
maximum THz energy yield of ^ 1.8 /iJ; 9^^ ^ — 7r/4 
leads to a minimum yield of 0.25 /iJ; the generic value 
^in = yields ^ 1 fj,3. THz spectra obtained for the two 



extremal cases are shown in Fig.[21Jc), cases (2) and (1), 
respectively. They agree remarkably well with the pre- 
dictions obtained from Eq. (j?]). 

We also observe a high sensitivity of the THz spec- 
tra towards detuning of the relative frequency ratio be- 
tween the two pump components [see case (3) and (4) in 
Fig. [2] . Such pulses with so-called incommensurate SH 
(frequency 2i'o is shifted by Sv) offer a straightforward 
control of the position of the maximum THz spectral 
intensity, without decreasing the total THz yield. Re- 
cently, this shift of the THz spectrum with Sv was pre- 
dicted [20| and demonstrated experimentally 21 1 . Again, 
we can use our linear interference framework to explain 
this effect. The responsible shifts of the ionization events 
can be computed approximatively (see Auxiliary mate- 
rial) as ojotn ss nTT — 2r(— 1)" sin {9 ± TmSv/vo). The SH 
frequency shift ±5^ makes the effective phase angle n 
dependent, and thus shifts the ionization events in time. 
Using our diffraction-grating analogy, the effect corre- 
sponds to increasing or decreasing the slit-distance along 
the grating. This frequency shift can be also calculated 
analytically under some further approximations (see the 
Auxiliary Material). 

Because the coefficients C„ are not sign definite [see 
Fig.[lla,b), the "secondary waves" are not in-phase], such 
non-equidistant grating can lead to a spatial shift of the 
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FIG. 3. THz generation by two-color pump pulses with in- 
commensurate frequencies, (a) THz spectrum vs. frequency 
shift Sv of the second harmonic computed from Eq. ([2| . The 
dashed white line shows the position of the spectral maxi- 
mum. (b),(c),(d) THz spectrum (color plot) as well as the 
maximum plasma density (white line) along the propagation 
axis according to 3D simulations for Sv — and ±20 THz, 
respectively. 



zeroth difFraction order, i.e., a frequency shift of the max- 
imum THz spectral intensity. As it is seen from Fig. HJc) 
and Fig. EJa) the THz frequency shift is not symmet- 
ric with respect to the sign of Siy. This effect occurs for 
pulses with finite duration (see Auxiliary Material). 

Again, our predictions from Eq. (|4]) are in good agree- 
ment with 3D numerical simulations. Incommensurate 
two-color pulses clearly produce the predicted THz spec- 
tral shift for negative SH detuning [see Fig. [Sjc)]. More- 
over, neither linear nor nonlinear propagation effects of 
the pump pulse due to focusing or laser-induced plasma 
seem to jeopardize this effect. The energy in the shifted 
THz pulse is ~ 1 /xJ, like in the corresponding commen- 
surate case [see Fig. El^b)]. On the other hand, Fig. ^d) 
reveals that, as expected, a positive SH dettming does 
not alter the THz spectral shape significantly, but the 
energy yield reduces by a factor of two to ~ 0.5 fiJ. 

Let us finally consider a more complex example of THz 
spectral engineering. Namely, we attempt to significantly 
increase the spectral width by using a field shape which 
optimizes the "amplitudes" C„ of just a few ionization 
events, suppressing the others. This leads to THz-to- 
mid-infrared supercontinuum, which is of great impor- 
tance for many applications. For this purpose, we resort 
to a three-color pump setup involving an optical para- 
metric amplifier (OPA). We assume a fundamental fre- 
quency vq = 375 THz (Ao = 800 nm) and signal and idler 
frequencies Vg and t'l obeying i^s + Vi = vq (t's = 0.55z^Oj 
z^i = 0.451^0). From Fig. IH^b), one can see that in such a 
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FIG. 4. Comparison of the THz characteristics produced 
by two-and three-color pump pulses, (a) Field shape E{t) 
(gray line), free electron velocity Vf{t) (blue line) as well as 
ionization amplitudes C„ (bars at ionization events) for the 
two-color pump (see Fig.[T}. (b) Same quantities for the three- 
color pump, e. g., for an OPA with frequencies vo = i^s + Vi 
{vs ~ 0.55!-'o, i^i ~ 0A5vo)- (c) THz energies U vs. propa- 
gation distance z according to 3D simulations and (d) cor- 
responding spatially integrated THz spectra at the position 
of the linear focus z — 5 mm for three-color (blue lines) and 
two-color (gray lines) cases. 



configuration only a few ionization events have large am- 
plitudes I Cn I ■ This immediately leads to a much broader 
THz spectrum, as it is obvious from the above discussed 
diffraction grating analogy. In Fig. IHc,d) the results of 
3D simulations for this case are shown, assuming that 
each of the three pump components contains one third of 
the total pulse energy of 300 fiJ while all other parame- 
ters are kept unchanged. 

Quite interestingly, the three-color configuration we 
use provides much larger values for the |C„| than the 
previous two-color scheme with comparable pump energy 
[cf. Fig.mja) and (b)]. In the LC limit, the resuUing THz 
yield is about 40 times larger. In our simulations, this in- 
creased THz yield is indeed visible in the beginning of the 
propagation z < 2.5 mm. However, upon further prop- 
agation saturation effects, mainly coming from plasma 
defocusing of the pump and depletion of neutral atoms, 
limit the total THz energy to a few /iJ. 



CONCLUSION AND OUTLOOK 

In summary, we have shown that THz emission in gases 
by muhi-color laser fields is intrinsically connected with 
the attosecond jumps of the plasma density induced by 
tunneling ionization. The THz spectrum is described by 
a superposition of the spectrally ultrabroad contributions 
from separate ionization events. 

This approach explains basic features of THz radiation 
such as the inverse dependence of the spectral width on 
the pulse duration and the shift of THz spectral max- 
ima for non-integer frequency offsets of two-color pumps. 
Moreover it allows to manipulate the THz waveform and 
to shape the THz spectrum by controlling the tunneling 
photo-ionization. In the present article we demonstrated 
only basic features of such a control, using frequencies 
and phases of two- or three-colour pump fields. We note 
that other degrees of freedom of light such as polarization 
or spatial shape of the pump pulses can also be used for 
THz control. Noteworthy, our approach is not limited to 
THz, far- and mid-infrared range but can be potentially 
extended to frequencies much higher than the pump one, 
i.e., to optical harmonics. 

We believe that this novel strategy to control fun- 
damental properties of THz pulses could become a ba- 
sic tool in optical technology: a ultra-broadband pulse 
shaper working in the range from THz to mid-infrared. 
Just as optical pulse shapers for femtosecond pulses are 
indispensable in ultrafast physics, the proposed THz con- 
trol can find promising applications in many fields of re- 
search, including, among others, THz time-domain spec- 
troscopy, chemical reaction dynamics, and optical signal 
processing. 



Methods 
Discrete ionization model 

The electron density p{t) is described by the equation 

dtP^WsTiE)[pat-pit)], (5) 

where pat denotes the neutral atomic density. 
We use a quasi-static tunneling ionization rate 
for hydrogen-like atoms [lO|, ll5|: Wst{E) = 
Aua{rH?-^[Ea/\E\] exp[-2(rj,)i-5i;j3|£;|], where 

Ea = m2g7(47reo)3ft4^ w^ = mq'^/{4:TTeo)'^h^, and 
rn = UAr/Uh- Ufi and UAr sue the ionization potentials 
of hydrogen and argon, respectively. 

Supposing that the ionization process occurs only near 
maxima of the pump electric field at times t„, one finds 
E{t) K Eo + E2{t - t„)2 with E2 - dttE{tn)/2 and 
Eq = E(tn). Then, assuming pat ^ p, we can evalu- 
ate WsTiE) w WsT[^(i„)]exp{-a2(t-t„)V4} where 



a^ = 8{rHy-^Ea\E2\mEo\^ and find 

p{t)=J2^PnH{t-tn), 



with 



-r2 (t— +„ a2 



(6) 



a{t - f„) 



H{t-tn)(x I e^' '"4*"' dr^W^O + erf 



Strictly speaking, a is different for every ionization event. 
The difference, however, is typically less than a few per- 
cent, so we assume here a "typical" value cr ^ 10 fs~^. 
The free electron current J is given by 



J{t)=q I [drPir)] Vfit)-e'<<^--%f{T) 



dr. 



which is equivalent to Eq. ([T|) (see [ll|, |15|, |22|). Using 
Eq. (ini), the LC J can be written as a sum over all contri- 
butions from separate ionization events J(i) = ^^ Jnit), 
and hence 



J{t) ^Y,q5p„H{t - t„) \^vf{t) 



piiu-t) 



"/ 



itn)]. (7) 



In the LC limit, assuming that the plasma source is con- 
tained in a small volume Ay, the field produced by 
the plasma current J{uj) in frequency domain is E'^ — 
gujj{uj) [27|, where g — i^^^ exp (iw^), c is the speed 
of light. So is the vacuum permittivity and r is the dis- 
tance between the plasma source and the point where the 
field is measured. 

Thus, in order to obtain Eq. ([2]), we have to multiply 
the Fourier transform of Eq. ([7]) by ujg. Further on, by 
performing the Fourier transform of H(t) exp (— 7^) and 
neglecting a slowly varying phase term exp {—2i^uj/a'^) 
we obtain Eq. ([3]). 



Comprehensive numerical model 

We use the unidirectional pulse propagation equation 
ISl |36J for linearly polarized pulses 



d.Eiuj) = z^B{u) -kl- klE{u^) + ^I^VMu) 



(8) 

Here, E(kx,ky, z,uj) is the relevant electric field compo- 
nent in Fourier domain with respect to the transverse co- 
ordinates {x,y) and time, k = ujn{u))/c is the wavenum- 
ber, c is the speed of light and n(w) is the linear re- 
fractive index of argon [37| . The nonlinear polarization 
'Pnl(i^) = ^Kcrr(w) -l- iJ{uj)/u} + «Jioss(w)/a; accounts 
for third-order nonlinear polarization PKcrrit), electron 
current J{t) and a loss term Jioss(0 due to photon ab- 
sorption during ionization. In our 3D numerical code. 



Eqs. dH), dl]), and ([5]) are solved using a standard spectral 
operator splitting scheme. Particular attention is paid to 
resolve all relevant time scales, from tens of attosecond 
(ionization steps) to a few picosecond (THz radiation). 
Although, strictly speaking, Eq. ([5|) does not describe the 
propagation of waves below the plasma frequency prop- 
erly, this plays only a minor role because the THz radi- 
ation is emitted from the leading plas ma front and thus 
does not propagate in the plasma [38|. Spectra obtained 
using Eq. (J8]) were found to be in good agreement with 
experimental results 
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Auxiliary Material 

We consider for simplicity the case of a two-color 
pulse with duration Tq having the square-shaped enve- 
lope E{t) = for \t\ > To/2 and 

E{t) = Ax cosjtJoi} + A2 cos{(2tJo + ^'^)i + (^} (9) 

for \t\ < Tq/2. Here Suj is a detuning from the second 
harmonic, which is assumed to be small compared to the 
main frequency {Suj <^ ujq). We also neglect the plasma 
decay 7 and assume r = A2/A1 <§:_ I. 

To determine the positions i„ of the tunneling ioniza- 
tion events (maxima of the electric field) along the time 
axis we equate the time derivative of Eq. ^ to zero, set 
tn = tno + Stn and expand the resulting expression up to 
first orders in Suj and Stn, yielding 

UJotn ~ mr — 2r(— 1)" [nnSuj/uJo cos0 + sin0) 

san7r-2r(-l)"sin(6' + n7r(5w/a;o), (10) 

where it is assumed that irnSoj/ujo <C 1. Using Eqs. 1^, 
(jlOp and neglecting again higher-order terms, we find 



Vf{tn) = — / E{T)dT w ai 



na2. 



where 



02 



: —3qAirsm6/2mujQ^ 
SgnAirSu! cosd/2mujQ. 



(11) 



(12) 
(13) 



We further assume that the amplitude Spn of each step 
in the modulation of the plasma density does not depend 
on the index n. For the square-shaped pulse envelope 
considered here this essentially means that saturation of 
the free electron density is neglected. Then, according to 
Eq. (4) of the main article 

E^{uj)oc{aiSi+a2S2), 
where 



(14) 



^1 



So 



^e-*", 



Y^ne'^*". 



(15) 
(16) 



Here, the summation runs over all ionization events. Fol- 
lowing Eq. (TTOl) . the exponent e*'^*" can be expressed 
as e-2*''""'"^+"'s+" for even n and e^'™'''"^-™^-'^ for 
odd n, where a; = uj/ujq is the normalized frequency and 
g+ — Tr{2r Suj cos 9 /ujQ + 1), g- — tt {2r Suj cos 9/ ujq — 1). 
Considering an odd number of ionization events N, being 
much larger than unity, we calculate the expressions in 
Eq. ([T4| by performing the summation from —M to M, 
where M = (A^ - l)/2. The term Si in Eq. ^ com- 
prises simple geometric series and gives after separate 
summation for odd and even n: 



Si 



„2^™ sin sin (Mg+Cd) ^ ^_2i™ sin 6 siu (Mg^uj) 
sin (g+uj) 



sin(g-w) 

(17) 
where M ^ 1, so that M -I- 1 w Af. A similar analy- 
sis is possible for the term 5*2 in Eq. (fT6|) . As before, 
one performs the sum for odd and even values of n sep- 
arately. We are then left with an arithmetico-geometric 
series [39j, i.e., for arbitrary q and integers mi, 1712: 



^ „ _ mi(g - 1) + g + g'»i+"'2+i [m2{q ^ 1) - 1] 

n—mi 



q'^^{q- 1)2 



After substituting Eq.([THl) into Eq. ([H]) we find 



(18) 



_i ( {M + 2) sin [Mg^uj] - M sin [(M + 2)g^uj] 
'"2\ sin^ {g_Gj) 

X exp {~2iruj sin 9)+ 
{M + 1) sin [{M - l)g+Cj] - (M - 1) sin [{M + l)g+Cj] 



sin^ {g+uj) 
X exp (2iruj sin 9) 



(19) 



Important features can be refound from the above ex- 
pressions. In particular, for Suj = we have 02 = and 
hence E'' (x Si. In addition, g+ = — g_ = tt and thus, 
by introducing St = tt/wq we obtain: 

f^ ■ r.^\sm (MujSt) 

S'iw2cos 2rsin6'— V ^ / . (20) 

\ UJQ J sm [ujSt] 

The equation given at the end of Sec. II in the main 
article is obtained from Eq. (f20l) by assuming 9 = Q and 
A^ » 1, that is, M ={N - l)/2 « N/2. 

For nonzero detuning Suj we can at least qualitatively 
reproduce the behavior presented in Fig. 2 (case 3) and 
in Fig. 3(a) of the main article, i.e., the shift of the max- 
ima of the THz spectra from zero frequency. The fre- 
quency shift obtained from Eq. ([TH), Eq. dTT]), Eq. p9)) 
is, however, symmetric with respect to the sign of (5a;, in 
contrast to the full LC model. The asymmetry reported 
in the main text [see, e.g.. Fig. 1(d)] disappears with 
increasing of the pulse duration and thus is induced by 
the presence of the Gaussian pulse envelope. 




5 OJO"/"0 




SoWt^O 




21) 0,25 0.3CJ<^/'^0 




o.y ' ' aj ^/'^o 



Typical form factors (a) \Si\^ , (b) 152 1^ and (c) the 



FIG. 5. _ 

field \E''\'^ depending on a) = uj/loo obtained from Eq. (I14|l . 
Eq. 1121), Eq. (Ull) for 6 = 7r/4, Af = 15, r = VO^, 5^M = 
0.0533 {5v = 20 THz). For comparison, in (d) the spectral 
shapes from Fig. 2(c) of the main article are repeated [blue 
curve corresponds to the case 3 (Su = — 20 THz) and red 
curve to the case 4 {5v = 20 THz)] . The insets in (c) and (d) 
show the corresponding pump pulse shapes. 
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